library(MASS)
library(haven)
library(Matrix)
library(lfe)
library(ggplot2)
library(RColorBrewer)
library(Hotelling)
library(did)

source("functions.R")

### run data_cleaning.R and classification.R first
load("data_classified.RData")

### K=2
# Table 5 of the manuscript
# compare treated and never-treated within type 1
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK2$gamma[temp$id]
temp <- temp[temp$type==1,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# compare treated and never-treated within type 2
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK2$gamma[temp$id]
temp <- temp[temp$type==2,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))
remove(temp)

# Table 6 of the manuscript
# compare type 1 and type 2
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK2$gamma[temp$id]
matrix(c(mean(temp$disd[temp$type==1]), sd(temp$disd[temp$type==1]),
         mean(temp$disd[temp$type==2]), sd(temp$disd[temp$type==2]),
         mean(temp$disd[temp$type==1])-mean(temp$disd[temp$type==2]),
         (var(temp$disd[temp$type==1])/length(temp$disd[temp$type==1])+
            var(temp$disd[temp$type==2])/length(temp$disd[temp$type==2]))^(0.5)),
       2,3)
matrix(c(mean(temp$ccbase[temp$type==1]), sd(temp$ccbase[temp$type==1]),
         mean(temp$ccbase[temp$type==2]), sd(temp$ccbase[temp$type==2]),
         mean(temp$ccbase[temp$type==1])-mean(temp$ccbase[temp$type==2]),
         (var(temp$ccbase[temp$type==1])/length(temp$ccbase[temp$type==1])+
            var(temp$ccbase[temp$type==2])/length(temp$ccbase[temp$type==2]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$type==1]), sd(temp$p_white[temp$type==1]),
         mean(temp$p_white[temp$type==2]), sd(temp$p_white[temp$type==2]),
         mean(temp$p_white[temp$type==1])-mean(temp$p_white[temp$type==2]),
         (var(temp$p_white[temp$type==1])/length(temp$p_white[temp$type==1])+
            var(temp$p_white[temp$type==2])/length(temp$p_white[temp$type==2]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$type==1]), sd(temp$p_hisp[temp$type==1]),
         mean(temp$p_hisp[temp$type==2]), sd(temp$p_hisp[temp$type==2]),
         mean(temp$p_hisp[temp$type==1])-mean(temp$p_hisp[temp$type==2]),
         (var(temp$p_hisp[temp$type==1])/length(temp$p_hisp[temp$type==1])+
            var(temp$p_hisp[temp$type==2])/length(temp$p_hisp[temp$type==2]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$type==1]), sd(temp$p_freelunch[temp$type==1]),
         mean(temp$p_freelunch[temp$type==2]), sd(temp$p_freelunch[temp$type==2]),
         mean(temp$p_freelunch[temp$type==1])-mean(temp$p_freelunch[temp$type==2]),
         (var(temp$p_freelunch[temp$type==1])/length(temp$p_freelunch[temp$type==1])+
            var(temp$p_freelunch[temp$type==2])/length(temp$p_freelunch[temp$type==2]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$type==1]), sd(temp$n_student[temp$type==1]),
         mean(temp$n_student[temp$type==2]), sd(temp$n_student[temp$type==2]),
         mean(temp$n_student[temp$type==1])-mean(temp$n_student[temp$type==2]),
         (var(temp$n_student[temp$type==1])/length(temp$n_student[temp$type==1])+
            var(temp$n_student[temp$type==2])/length(temp$n_student[temp$type==2]))^(0.5)),
       2,3)
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,2)))

### K=3
# Table 5 of SA
# compare treated and never-treated within type 1
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK3$gamma[temp$id]
temp <- temp[temp$type==1,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# compare treated and never-treated within type 2
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK3$gamma[temp$id]
temp <- temp[temp$type==2,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# Table 6 of SA
# compare treated and never-treated within type 3
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK3$gamma[temp$id]
temp <- temp[temp$type==3,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# Table 9 of SA
# compare type 1, 2 and 3
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK3$gamma[temp$id]
matrix(c(mean(temp$disd[temp$type==1]), sd(temp$disd[temp$type==1]),
         mean(temp$disd[temp$type==2]), sd(temp$disd[temp$type==2]),
         mean(temp$disd[temp$type==3]), sd(temp$disd[temp$type==3]),
         mean(temp$disd[temp$type==1])-mean(temp$disd[temp$type==2]),
         (var(temp$disd[temp$type==1])/length(temp$disd[temp$type==1])+
            var(temp$disd[temp$type==2])/length(temp$disd[temp$type==2]))^(0.5)),
       2,4)
matrix(c(mean(temp$ccbase[temp$type==1]), sd(temp$ccbase[temp$type==1]),
         mean(temp$ccbase[temp$type==2]), sd(temp$ccbase[temp$type==2]),
         mean(temp$ccbase[temp$type==3]), sd(temp$ccbase[temp$type==3]),
         mean(temp$ccbase[temp$type==1])-mean(temp$ccbase[temp$type==2]),
         (var(temp$ccbase[temp$type==1])/length(temp$ccbase[temp$type==1])+
            var(temp$ccbase[temp$type==2])/length(temp$ccbase[temp$type==2]))^(0.5)),
       2,4)
matrix(c(mean(temp$p_white[temp$type==1]), sd(temp$p_white[temp$type==1]),
         mean(temp$p_white[temp$type==2]), sd(temp$p_white[temp$type==2]),
         mean(temp$p_white[temp$type==3]), sd(temp$p_white[temp$type==3]),
         mean(temp$p_white[temp$type==1])-mean(temp$p_white[temp$type==2]),
         (var(temp$p_white[temp$type==1])/length(temp$p_white[temp$type==1])+
            var(temp$p_white[temp$type==2])/length(temp$p_white[temp$type==2]))^(0.5)),
       2,4)
matrix(c(mean(temp$p_hisp[temp$type==1]), sd(temp$p_hisp[temp$type==1]),
         mean(temp$p_hisp[temp$type==2]), sd(temp$p_hisp[temp$type==2]),
         mean(temp$p_hisp[temp$type==3]), sd(temp$p_hisp[temp$type==3]),
         mean(temp$p_hisp[temp$type==1])-mean(temp$p_hisp[temp$type==2]),
         (var(temp$p_hisp[temp$type==1])/length(temp$p_hisp[temp$type==1])+
            var(temp$p_hisp[temp$type==2])/length(temp$p_hisp[temp$type==2]))^(0.5)),
       2,4)
matrix(c(mean(temp$p_freelunch[temp$type==1]), sd(temp$p_freelunch[temp$type==1]),
         mean(temp$p_freelunch[temp$type==2]), sd(temp$p_freelunch[temp$type==2]),
         mean(temp$p_freelunch[temp$type==3]), sd(temp$p_freelunch[temp$type==3]),
         mean(temp$p_freelunch[temp$type==1])-mean(temp$p_freelunch[temp$type==2]),
         (var(temp$p_freelunch[temp$type==1])/length(temp$p_freelunch[temp$type==1])+
            var(temp$p_freelunch[temp$type==2])/length(temp$p_freelunch[temp$type==2]))^(0.5)),
       2,4)
matrix(c(mean(temp$n_student[temp$type==1]), sd(temp$n_student[temp$type==1]),
         mean(temp$n_student[temp$type==2]), sd(temp$n_student[temp$type==2]),
         mean(temp$n_student[temp$type==3]), sd(temp$n_student[temp$type==3]),
         mean(temp$n_student[temp$type==1])-mean(temp$n_student[temp$type==2]),
         (var(temp$n_student[temp$type==1])/length(temp$n_student[temp$type==1])+
            var(temp$n_student[temp$type==2])/length(temp$n_student[temp$type==2]))^(0.5)),
       2,4)
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,2)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,3)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(2,3)))


### K=4
# Table 7 of SA
# compare treated and never-treated within type 1
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK4$gamma[temp$id]
temp <- temp[temp$type==1,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# compare treated and never-treated within type 2
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK4$gamma[temp$id]
temp <- temp[temp$type==2,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# Table 8 of SA
# compare treated and never-treated within type 3
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK4$gamma[temp$id]
temp <- temp[temp$type==3,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)
print(hotelling.test(ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ treated,
                     data=temp))

# compare treated and never-treated within type 4
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK4$gamma[temp$id]
temp <- temp[temp$type==4,]
temp$treated <- as.factor(temp$E <= T)
matrix(c(mean(temp$ccbase[temp$treated==TRUE]), sd(temp$ccbase[temp$treated==TRUE]),
         mean(temp$ccbase[temp$treated==FALSE]), sd(temp$ccbase[temp$treated==FALSE]),
         mean(temp$ccbase[temp$treated==TRUE])-mean(temp$ccbase[temp$treated==FALSE]),
         (var(temp$ccbase[temp$treated==TRUE])/length(temp$ccbase[temp$treated==TRUE])+
            var(temp$ccbase[temp$treated==FALSE])/length(temp$ccbase[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_white[temp$treated==TRUE]), sd(temp$p_white[temp$treated==TRUE]),
         mean(temp$p_white[temp$treated==FALSE]), sd(temp$p_white[temp$treated==FALSE]),
         mean(temp$p_white[temp$treated==TRUE])-mean(temp$p_white[temp$treated==FALSE]),
         (var(temp$p_white[temp$treated==TRUE])/length(temp$p_white[temp$treated==TRUE])+
            var(temp$p_white[temp$treated==FALSE])/length(temp$p_white[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_hisp[temp$treated==TRUE]), sd(temp$p_hisp[temp$treated==TRUE]),
         mean(temp$p_hisp[temp$treated==FALSE]), sd(temp$p_hisp[temp$treated==FALSE]),
         mean(temp$p_hisp[temp$treated==TRUE])-mean(temp$p_hisp[temp$treated==FALSE]),
         (var(temp$p_hisp[temp$treated==TRUE])/length(temp$p_hisp[temp$treated==TRUE])+
            var(temp$p_hisp[temp$treated==FALSE])/length(temp$p_hisp[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$p_freelunch[temp$treated==TRUE]), sd(temp$p_freelunch[temp$treated==TRUE]),
         mean(temp$p_freelunch[temp$treated==FALSE]), sd(temp$p_freelunch[temp$treated==FALSE]),
         mean(temp$p_freelunch[temp$treated==TRUE])-mean(temp$p_freelunch[temp$treated==FALSE]),
         (var(temp$p_freelunch[temp$treated==TRUE])/length(temp$p_freelunch[temp$treated==TRUE])+
            var(temp$p_freelunch[temp$treated==FALSE])/length(temp$p_freelunch[temp$treated==FALSE]))^(0.5)),
       2,3)
matrix(c(mean(temp$n_student[temp$treated==TRUE]), sd(temp$n_student[temp$treated==TRUE]),
         mean(temp$n_student[temp$treated==FALSE]), sd(temp$n_student[temp$treated==FALSE]),
         mean(temp$n_student[temp$treated==TRUE])-mean(temp$n_student[temp$treated==FALSE]),
         (var(temp$n_student[temp$treated==TRUE])/length(temp$n_student[temp$treated==TRUE])+
            var(temp$n_student[temp$treated==FALSE])/length(temp$n_student[temp$treated==FALSE]))^(0.5)),
       2,3)

# Table 10 of SA
# compare type 1,2,3 and 4
temp <- LutzData[KmeansData$time==1,]
temp$type <- kmeansK4$gamma[temp$id]
matrix(c(mean(temp$disd[temp$type==1]), sd(temp$disd[temp$type==1]),
         mean(temp$disd[temp$type==2]), sd(temp$disd[temp$type==2]),
         mean(temp$disd[temp$type==3]), sd(temp$disd[temp$type==3]),
         mean(temp$disd[temp$type==4]), sd(temp$disd[temp$type==4])),
       2,4)
matrix(c(mean(temp$ccbase[temp$type==1]), sd(temp$ccbase[temp$type==1]),
         mean(temp$ccbase[temp$type==2]), sd(temp$ccbase[temp$type==2]),
         mean(temp$ccbase[temp$type==3]), sd(temp$ccbase[temp$type==3]),
         mean(temp$ccbase[temp$type==4]), sd(temp$ccbase[temp$type==4])),
       2,4)
matrix(c(mean(temp$p_white[temp$type==1]), sd(temp$p_white[temp$type==1]),
         mean(temp$p_white[temp$type==2]), sd(temp$p_white[temp$type==2]),
         mean(temp$p_white[temp$type==3]), sd(temp$p_white[temp$type==3]),
         mean(temp$p_white[temp$type==4]), sd(temp$p_white[temp$type==4])),
       2,4)
matrix(c(mean(temp$p_hisp[temp$type==1]), sd(temp$p_hisp[temp$type==1]),
         mean(temp$p_hisp[temp$type==2]), sd(temp$p_hisp[temp$type==2]),
         mean(temp$p_hisp[temp$type==3]), sd(temp$p_hisp[temp$type==3]),
         mean(temp$p_hisp[temp$type==4]), sd(temp$p_hisp[temp$type==4])),
       2,4)
matrix(c(mean(temp$p_freelunch[temp$type==1]), sd(temp$p_freelunch[temp$type==1]),
         mean(temp$p_freelunch[temp$type==2]), sd(temp$p_freelunch[temp$type==2]),
         mean(temp$p_freelunch[temp$type==3]), sd(temp$p_freelunch[temp$type==3]),
         mean(temp$p_freelunch[temp$type==4]), sd(temp$p_freelunch[temp$type==4])),
       2,4)
matrix(c(mean(temp$n_student[temp$type==1]), sd(temp$n_student[temp$type==1]),
         mean(temp$n_student[temp$type==2]), sd(temp$n_student[temp$type==2]),
         mean(temp$n_student[temp$type==3]), sd(temp$n_student[temp$type==3]),
         mean(temp$n_student[temp$type==4]), sd(temp$n_student[temp$type==4])),
       2,4)
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,2)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,3)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(1,4)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(2,3)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(2,4)))
print(hotelling.test(disd + ccbase + p_white + p_hisp + p_freelunch + n_student
                     ~ type,
                     data=temp,pair=c(3,4)))